﻿using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;

namespace Pyramid
{
    public class GImage
    {
        public int Width;
        public int Height;
        public int Depth;
        public double[, ,] data;
        public GImage(int w, int h)
        {
            //create an empty gimage
            Width = w;
            Height = h;
            Depth = 3;
            data = new double[Width, Height, Depth];
        }
  
        public GImage(GImage newimage)
        {
            //create a gimage from another gimage
            Width = newimage.Width;
            Height = newimage.Height;
            Depth = newimage.Depth;
            data = new double[Width, Height, Depth];
            data = newimage.data;
        }
        public GImage(Bitmap newimage)
        {
            //create a gimage from a bitmap image
            Width = newimage.Width;
            Height = newimage.Height;
            Depth = 3;
            data = new double[Width, Height, Depth];
            for(int i=0;i<Width;i++)
                for (int j = 0; j < Height; j++)
                {
                    data[i, j, 0] = newimage.GetPixel(i, j).R;
                    data[i, j, 1] = newimage.GetPixel(i, j).G;
                    data[i, j, 2] = newimage.GetPixel(i, j).B;
                }
            
        }
        public static GImage operator +(GImage image1, GImage image2)
        {
            //only valid if the two images are of the same size
            if (image1.Width == image2.Width && image1.Height == image2.Height)
            {
                GImage result = new GImage(image1.Width, image1.Height);
                for (int i = 0; i < result.Width; i++)
                    for (int j = 0; j < result.Height; j++)
                    {
                        double[] r = new double[result.Depth];
                        for (int k = 0; k < result.Depth; k++)
                            r[k] = image1.data[i, j, k] + image2.data[i, j, k];
                        result.SetPixel(i, j, r);
                    }
                return result;
            }
            else
                return null;
        }
        public static GImage operator -(GImage image1, GImage image2)
        {
            //only valid if the two images are of the same size
            if (image1.Width == image2.Width && image1.Height == image2.Height)
            {
                GImage result = new GImage(image1.Width, image1.Height);
                for (int i = 0; i < result.Width; i++)
                    for (int j = 0; j < result.Height; j++)
                    {
                        double[] r = new double[result.Depth];
                        for (int k = 0; k < result.Depth; k++)
                            r[k] = image1.data[i, j, k] - image2.data[i, j, k];
                        result.SetPixel(i, j, r);
                    }
                return result;
            }
            else
                return null;
        }
        public static GImage operator *(double s, GImage image1)
        {
            //scalar multiplication
            GImage result = new GImage(image1.Width, image1.Height);
            for (int i = 0; i < result.Width; i++)
                for (int j = 0; j < result.Height; j++)
                {
                    double[] r = new double[result.Depth];
                    for (int k = 0; k < result.Depth; k++)
                        r[k] = s * image1.data[i, j, k];
                    result.SetPixel(i, j, r);
                }
            return result;
        }
        public bool SetPixel(int row, int col, double[] val)
        {
            //set the pixel at location (row,col) to val
            if (row < Width && col < Height && val.Length==Depth)
            {
                for (int i = 0; i < Depth; i++)
                    data[row, col, i] = val[i];
                return true;
            }
            else
                return false;
        }
        public GImage Normalize()
        {
            //change the range of the values from whatever it is to [0,255]
            //min<x<max ==> 0<255*(x-min)/(max-min)<255
            GImage result = new GImage(Width, Height);
            double[,] m = MinMax(1);

            
            double[] r = new double[Depth];
            int i, j, k;
            

            

            for (i = 0; i < Width; i++)
                for (j = 0; j < Height; j++)
                {
                    for (k = 0; k < Depth; k++)
                    {
                        r[k] = 255.0*(data[i,j,k] - m[0, k]) / (m[1, k] - m[0, k]);
                    }
                    result.SetPixel(i, j, r);
                }
           
            return result;
        }
        public double[,] MinMax(int mode)
        {
            //return the mazimum and minimum of the image in a 2D array
            //the first row is min and the second row is max
            //mode=0 find different min and max for different components
            //mode=1 find one min and max for all the components
            double[,] m=new double[2,Depth];
            double min, max;
            int i,j,k;
            for (i = 0; i < Depth; i++)
                m[0, i] = 100000000;//big number
            min = 100000000;
            for (i = 0; i < Depth; i++)
                m[1, i] = -100000000;//small number
            max = -10000000;
            for(i=0;i<Width;i++)
                for (j = 0; j < Height; j++)                
                    for (k = 0; k < Depth; k++)
                    {
                        if (data[i, j, k] < m[0, k])
                            m[0, k] = data[i, j, k];
                        if (data[i, j, k] > m[1, k])
                            m[1, k] = data[i, j, k];

                        if (data[i, j, k] < min)
                            min = data[i, j, k];
                        if (data[i, j, k] > max)
                            max = data[i, j,k];

                    }
            if (mode == 1)
            {
                for (k = 0; k < Depth; k++)
                {
                    m[0, k] = min;
                    m[1, k] = max;
                }
            }
            return m;
        }
        public List<double[]> FindNeighborhood(int[] center, int NeighborhoodSize)
        {
            int i,j,k;

            List<double[]> n = new List<double[]>();//the pixel values in the neighborhood
            // the new pixel and its color
            double[] pixel;

            if ((center[0]) < 0 || (center[0]) >= Width || (center[1]) < 0 || (center[1]) >= Height)
                return n;
            //the first value is always the center
            pixel = new double[Depth];
            for (i = 0; i < pixel.Length; i++)
                pixel[i] = data[center[0], center[1], i];

            n.Add(pixel);
            //find the neighbors
            for (i = -NeighborhoodSize / 2; i <= NeighborhoodSize / 2; i++)
            {
                
                if ((center[0] + i) >= 0 && (center[0] + i) <Width)
                {
                    for (j = -NeighborhoodSize / 2; j <= NeighborhoodSize / 2; j++)
                    {
                        if ((Math.Abs(i)+Math.Abs(j))!=0 && (center[1] + j) >= 0 && (center[1] + j) < Height)
                        {
                            pixel = new double[3];
                            for (k = 0; k < pixel.Length; k++)
                                pixel[k] = data[center[0] + i, center[1] + j, k];
                            n.Add(pixel);
                        }
                    }
                }
            }

            return n;
        }
        public double StructuralSimilarity(GImage image2, int NeighborhoodSize, int offsetx=0,int offsety=0)
        {
            double s;
            double[] SSIM = new double[image2.Depth];
            int i, j;
            int[] center = new int[2];
            double L = 255;
            Statistics MySt = new Statistics();
            //make sure ure not going out of the image
            if (NeighborhoodSize > image2.Width)
                NeighborhoodSize = image2.Width;
            if (NeighborhoodSize > image2.Height)
                NeighborhoodSize = image2.Height;

            for (i = 0; i < SSIM.Length; i++)
                SSIM[i] = 0;
            //march both images pixel by pixel
            for (i = 0; i + offsetx < image2.Width; i++)
                for (j = 0; j + offsety < image2.Height; j++)
                {
                    List<double[]> Neighborhood1 = new List<double[]>();
                    List<double[]> Neighborhood2 = new List<double[]>();
                    center[0] = i; center[1] = j;
                    Neighborhood1 = FindNeighborhood( center, NeighborhoodSize);
                    center[0] = i+offsetx; center[1] = j+offsety;
                    Neighborhood2 = image2.FindNeighborhood(center, NeighborhoodSize);
            
                    //1- build gaussian models for fine image neighborhood
                    double[] ms = MySt.FindMean(Neighborhood1);
                    double[] vs = MySt.FindVariance(Neighborhood1);

                    //2- build gaussian models for coarse image neighborhood
                    double[] mc = MySt.FindMean(Neighborhood2);
                    double[] vc = MySt.FindVariance(Neighborhood2);


                    double[] cov=MySt.FindCovariance(Neighborhood1,Neighborhood2);
                    //3-find the similarity
                  
                    
                    double k1 = 0.0001;//small constant <<1
                    double c1 = Math.Pow(k1*L, 2);
                  
                    
                    double k2 = 0.0001;//small constant <<1
                    double c2 = Math.Pow(k2 * L,2);

                    

                    for (int k = 0; k < SSIM.Length; k++)
                        SSIM[k] += ((2 * ms[k] * mc[k] + c1)*(2*cov[k]+c2)) /( (Math.Pow(ms[k], 2) + Math.Pow(mc[k], 2) + c1)* (Math.Pow(vs[k], 2) + Math.Pow(vc[k], 2) + c2));

                    
                }
            for (i = 0; i < SSIM.Length; i++)
                SSIM[i] = SSIM[i]/(image2.Width*image2.Height);
            s = 0;
            for (i = 0; i < SSIM.Length; i++)
                s+=SSIM[i];

            return s/SSIM.Length;
        }
        public Point Register(GImage image2,int SearchWindowSize)
        {
            //search for the subset of image2 which matches the best with this image
           
            Point r = new Point();//location of the best match
            int i, j,nsize;
            double max_s,curr_s;//the maximum similarity and the current similarity
            max_s = -1;
            //make sure ure not going out of the image
            if (SearchWindowSize > image2.Width)
                SearchWindowSize = image2.Width;
            if (SearchWindowSize > image2.Height)
                SearchWindowSize = image2.Height;
            nsize = 5;
            for(i=0;i<SearchWindowSize;i++)
                for (j = 0; j < SearchWindowSize; j++)
                {
                    //find the similarity between this image and the subset of image2 that starts from the location (i,j)
                    curr_s = StructuralSimilarity(image2, nsize,i,j);
                    if (curr_s > max_s)
                    {
                        max_s = curr_s;
                        r.X = i;
                        r.Y = j;
                    }
                }
            return r;
        }
        public Bitmap ToBitmap()
        {
            //convert the image to bitmap for display purposes
            Bitmap image = new Bitmap(Width, Height);
            int[] r = new int[3];

            for(int i=0;i<Width;i++)
                for (int j = 0; j < Height; j++)
                {
                    for (int k = 0; k < 3; k++)
                    {
                        r[k] = (int)(data[i, j, k]);
                        //push to the valid range of rgb: [0,255]
                        r[k] = r[k] < 255 ? r[k] : 255;
                        r[k] = r[k] > 0 ? r[k] : 0;
                    }
                    image.SetPixel(i,j, Color.FromArgb(r[0], r[1], r[2]));
                }
            return image;
        }
    }
}
